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The ABC-stacked A^-layer-graphene family of two-dimensional electron systems is described at low 
energies by two remarkably flat bands with Bloch states that have strongly momentum-dependent 
phase differences between carbon 7r-orbital amplitudes on different layers, and large associated mo- 
mentum space Berry phases. These properties are most easily understood using a simplified model 
with only nearest-neighbor inter-layer hopping which leads to gapless semiconductor electronic struc- 
ture, with dispersion in both conduction and valence bands. We report on a study of the elec- 
tronic band structures of trilayers which uses ah initio density functional theory and fc • p theory to 
fit the parameters of a 7r-band tight-binding model. We find that when remote interlayer hopping is 
retained, the triple Dirac point of the simplified model is split into three single Dirac points located 
along the three KM directions. External potential differences between top and bottom layers are 
strongly screened by charge transfer within the trilayer, but still open an energy gap at overall 
neutrality. 



PACS numbers: 73.43.Cd, 71.15.-m, 71.20.-b, 81.05.ue 

I. INTRODUCTION 

Success^ in isolating nearly perfect monolayer and few 
layer sheets from bulk graphite, along with progress^ in 
the epitaxial growth of few-layer graphene samples, has 
led to an explosion of experimental and theoretical^ri in- 
terest in this interesting class of quasi-two-dimensional 
electron systems (2DES's). Unique aspects of the elec- 
tronic structure of graphene based 2DES's have raised 
a number of new fundamental physics issues and raised 
hope for applications. 

Monolayer graphene has a honeycomb lattice structure 
and is a gapless semiconductor. Hopping between its 
equivalent A and B sublatticcs gives rise to a masslcss 
Dirac fermion band structure with J = 1 chirality when 
the sublattice degree of freedom is treated as a pseu- 
dospin. In this paper we will find it useful to view the 
quantum two-level degree of freedom associated with two 
sublattice sites as a pseudospin in the multi-layer case as 
well. In AB-stacked graphene bilayers, for example, elec- 
trons on the A2 and Bi sublatticcs arc repelled from the 
Fermi level by a direct interlayer tunneling process with 
energy 71, leaving^ only states that are concentrated on 
the Ai and B2 sites in the low-energy band-structure pro- 
jection. When direct hopping between Ai and B2 sites is 
neglected, the two-step hopping process via high energy 
sites leads to conduction and valence band dispersions 
and to a pseudospin chirality that is doubled, i.e. to 
a phase difference between sublattice projections which 
is proportional to 2(/)p where (j)p is the two-dimensional 
momentum orientation. Pseudospin chirality has a sub- 
stantial influence on interaction physics^ in both single- 
layer and bilayer graphene, and through the associated 
momentum space Berry phases also on Landau quantiza- 
tion and the integer quantum Hall effect 4ii2rJ^. Because 



the two low-energy sublatticcs in bilayer graphene are lo- 
cated on opposite layers it is possible to introducei^"— a 
gapiSi"— in the electronic structure simply by using gates 
to induce a difference in electric potential between lay- 
ers. According to some theories a small gap could even 
emerge spontancouslj«2i"— in neutral graphene bilayers 
with weak disorder because of layer inversion symmetry- 
breaking. 

Graphene bilayer 2DES's are quite distinct from sin- 
gle layer 2DES's because of their flatter band dispersion 
and the possibility of using external potentials to create 
gaps. Among all stacking possibilities, only the ABC ar- 
rangement (see below) maintains the following features 
that make Bernal bilayer electronic structure interesting 
in thicker iV-layer films: (i) there are two low-energy sub- 
lattice sites, implying that a two band model provides a 
useful tool to describe its physics; (ii) the low energy 
sublattice sites are localized in the outermost layers, at 
Ai and B^r, and can be separated energetically by an 
electric field perpendicular to the film; (iii) hopping be- 
tween low energy sites via high energy states is an iV- 
stcp process which leads to dispersion in conduction 
and valence bands, sublattice pseudospin chirality N and 
Berry phase A^tt. The low-energy bands are increasingly 
flat for larger A'^, at least when weak remote hopping 
processes are neglected, and the opportunity for interest- 
ing interaction and disorder physics is therefore stronger. 
Consequently, in the simplified chiral model, the density- 
of-states D{E) ^ i<;(2-JV)/Af (diverges as E approaches 
zero for N > 2 whereas it remains finite for N = 2 and 
vanishes for A'^ = 1. These properties also have some rel- 
evance to more general stacking arrangements since the 
low energy Hamiltonian of a multilayer with any type of 
stacking can always be chiral-decomposed to a direct sum 
of ABC-stacked layers.— 

ABC-stacked multilayers are the chiral generalizations 
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of monolayer and Bernal bilayer graphene, and we re- 
fer them collectively as the chiral 2D electron system 
(C2DES) family. We believe that they are likely to prove 
to be fertile ground for new physics. As a first step in the 
exploration of these materials we report in this paper on 
an effort to characterize the way in which the chirality N 
bands of an A^-layer C2DES arc altered by remote hop- 
ping processes neglected in the simplified model, focus- 
ing on the iV = 3 trilayer case. We use ab initio density 
functional theory (DFT) calculations, combined with a 
k ■ p expansion of the low-energy bands near the Dirac 
point, to fit the parameters of a phenomenological tight- 
binding method (PTBM) for the 7r-bands of multilayer 
graphene ji^i^^"— We find that details of the low-energy 
band dispersion can be used to fix rather definite values 
for the model's remote intcr-laycr hopping parameters. 

Our paper is organized as follows. In section II we 
first sketch the derivation of the low energy effective 
band Hamiltonian of a trilayer, reserving details to an 
Appendix and explain how the interlayer hopping pa- 
rameters infiuence the shape of constant energy surfaces. 
The values for these parameters obtained by fitting to 
our DFT calculations are surprisingly different from the 
values for the analogous hopping parameters in Bernal 
stacked layers, and are not yet available from experi- 
ment. In Section II we also discuss the evolution of con- 
stant energy surface pockets with energy, concentrating 
on the Lifshitz transitions at which pockets combine, in 
terms of Berry phase considerations and a competition 
between chiral dispersion and trigonal warping. In sec- 
tion III we use DFT to estimate the dependence of the 
trilayer energy gap on the external potential difference 
between top and bottom layers and compare with pre- 
dictions based on the simplified two-band model. The 
simplified model picture is readily extended to higher N 
and we use it to discuss trends in thicker ABC multilay- 
ers. Finally, we conclude in Section IV with a discussion 
of how Berry phases modify the integer quantum Hall 
effect and weak localization in C2DES's and with some 
speculations on the role of electron-electron interactions 
in these two-dimensional materials. 



II. EFFECTIVE MODEL AND BAND 
STRUCTURE 

A. Low Energy Effective Model 

In ABC-stacked graphene layers, each layer has in- 
equivalent triangular A and B sublatticcs. As illustrated 
in Fig.[TJa), each adjacent layer pair forms a AB-stackcd 
bilayer with the upper B sublattice directly on top of the 
lower A sublattice, and the upper A above the center of a 
hexagonal plaquette of the layer below. Our microscopic 
analysis uses the categorization of interlayer hopping pro- 
cesses illustrated in Fig. [IJb), which is analogous to the 
Slonczewski-Wciss-McClure (SWM) paramctrization of 
the tight binding model of bulk graphite with the Bernal 



stacking order Following convention 70 and 71 describe 
nearest neighbor intralayer and interlayer hopping re- 
spectively, 73 represents hopping between the low energy 
sites of a AB-stacked bilayer {i.e. Ai (z = 1,2)), 

74 couples low and high energy sites located on different 
layers {i.e. Ai <-> Aj+i and Bi o Bi+i {i = 1,2)). We 
use 72 to denote direct hopping between the trilayer low 
energy sites, and S as the on-site energy difference of Ai 
and i?3 with respect to the high energy sites. 75 and 75 
correspond to the presumably weaker couplings Bi A3 
and Si -!->■ S3 (S = A,B), respectively and Ui is used to 
denote the average potential of the ith layer. 
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FIG. 1: (Color online) (a) Lattice structure of ABC-stacked 
graphene trilayer; blue/cyan/green indicate links on the 
top/middle/bottom layers while purple/red distinguish the 
A/B sublattices. (b) Schematic of the unit cell of ABC- 
stacked graphene trilayer and the most important interlayer 
hopping processes. 

The massless Dirac- Weyl quasiparticles of monolayer 
graphene are described by a fe • p Hamiltonian, 
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TT 
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where tt — + ipy and ^ = +( — ) for valley K{K'). 
(In the rest of the paper we focus on bands near Bril- 
louin zone corner K] the general result can be obtained 
by setting Px to Cpx-) The trilayer 7r-bands are the di- 
rect produce of three sets of monolayer bands, modi- 
fied by the various interlayer coupling processes identified 
above. In a representation of sublattice sites in the or- 
der Ai, B^^Bi, A2, B2, A3, the trilayer Hamiltonian near 
valley K can then be expressed in the form: 
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where 



y/?>ayi/2h and a = 0.246nm. 
The identification of Ai and B3 as the low-energy sub- 
lattice sites is made by neglecting the weaker remote in- 
terlayer hopping processes and setting tt — 0. We treat 
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coupling between the low and high-energy subspaces per- 
turbatively by writing the trilayer Greens function as 
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where the indices 1 and 2 denote the 2x2 low-energy 
block and the 4x4 high-energy block respectively. We 
then solve the Schrodinger equation, ij}\ovi = 0, by 

using the block matrix inversion rule {A~^)ii = [An — 
Ai2{A22)~^ A2i)~^ to obtain 



((iJll - e) - Hl2(i?22 - e)"'i?2l)V'low(A„B3) =0. (4) 

Since we are interested in the low-energy part of the spec- 
trum we can view e as small compared to i?22- Expanding 
Eq. (|4]) to first order in e, we find that (iJoff — e)'*/'iow = 0, 
where 

H,s = (1 + Hi2{H22r^H2iy\Hii - HMH22)-'H2i) .(5) 

The terms in the second parenthesis capture the lead- 
ing hopping processes between low energy sites, includ- 
ing virtual hopping via high-energy states, while the first 
parenthesis captures an energy scale renormalization by 
a factor of order 1 — due to higher-order pro- 

cesses which we drop except in the terms which arise from 
an external potential. 

Using Eq. ([5]) we find that for ABC trilayer graphene 
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Here we have chosen tan(/3p = Py/px, defined Ud = 
(ui — uy,)/2 and Ua = {ui + u^)/2 — W2, and neglected 
an overall energy scale associated with the external po- 
tentials, do is the identity matrix and the di's are Pauli 
matrices acting on the low-energy pseudospin. We have 
retained leading terms with cubic, quadratic, and con- 
stant dispersions, which are due respectively to three- 
step, two-step, and one-step hopping processes between 
low energy sites. For trilayer graphene, the linear term 
is absent because the one step hopping (72) is normal 
to the 2D space and therefore independent of momen- 
tum. Hc\x is the only term which appears in the effective 



Hamiltonian in the simplified model with only nearest 
neighbor inter-layer tunneling. This term has pseudospin 
chirality J = 3 and dominates at larger values of p. It 
reflects coupling between low energy sites via a sequence 
of three nearest neighbor intralayer and interlayer hop- 
ping events, i/tr is proportional to ax and, because it is 
isotropic in 2D momentum space, is responsible for trig- 
onal warping of constant energy surfaces when combined 
with the J = 3 chiral term. Notice that the direct hop- 
ping 72 process opens a small gap at the K points so 
that -fftr vanishes at finite p if 72 is positive. H^ arises 
from a weaker coupling between low energy and high en- 
ergy states that is present in bilayers and for any > 1 
multilayer system. This term in the effective Hamilto- 
nian preserves layer inversion symmetry, -ffgap captures 
the external potential processes which break layer inver- 
sion symmetry and introduce a gap between electron and 
hole bands. The possibility of opening a gap with an ex- 
ternal potential is unique to ABC stacked multilayers, 
increasing the possibility that they could be useful ma- 
terials for future semiconductor devices. The strength of 
the gap term decreases with increasing momentum (since 
vqp <C 71 ) so that the gap around K has a Mexican hat 
shape, as we will discuss later. H'^ is non-zero when the 
potential of the middle layer deviates from the average 
of the potentials on the outermost layers. Unlike i?gap, 
this term preserves the layer inversion symmetry and is 
not responsible for an energy gap. A non-zero H'^ is rel- 
evant when the electric fields in the two inter-layer re- 
gions are different. Further discussion on the derivation 
of this effective Hamiltonian and on the physical mean- 
ing of the various terms can be found in the Appendix. 
Note that for strict consistency the constant terms 5 and 
72/2 should be accompanied by the factor 1 — (wop)^/7i 
based on Eq. which does appear in i^gap- However we 
ignore this factor because 5 and 72/2 are already small. 



B. AB Initio Density Functional Theory 
Calculations 

We have performed ah initio DFT calculations^ for 
an isolated graphene trilayer in the absence of a trans- 
verse external electric field which induces an electric 
potential difference between the layers. (DFT calcula- 
tions in the presence of electric fields will be discussed in 
the next section.) Our electronic structure calculations 
were performed with plane wave basis sets and ultra- 
soft pseudopotentials^. The local density approxima- 
tion (LDA) was used for the exchange and correlation 
potential. We fixed the interlayer separation to 0.335 
nm and placed bulk trilayer graphene in a supercell with 
a 40 nm vacuum region, large enough to avoid intercell 
interactions. A21x21xl k-point mesh in the full 
supercell Brillouin zone (FBZ) was used with a 408 eV 
kinetic energy cut-off. The calculations were tested for 
large k-point meshes in the FBZ and large energy cut- 
offs for convergence studies. Fig.[2]shows the DFT energy 
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FIG. 2: (Color online) Band structure of ABC-stacked 
graphene trilayers in the absence of an external electric field. 
The zero of energy in this plot is at the Fermi energy of a 
neutral trilayer. Notice the single low-energy band with ex- 
tremely flat dispersion near the K point. 

band structure of ABC stacked trilayer graphene in the 
absence of an external electric field. The low energy band 
dispersion is nearly cubic at the two inequivalent corners 
K and K' of the hexagonal Brillouin zone, as predicted 
by the 7r-orbital tight-binding and continuum model phe- 
nomenologies. The conduction and valence bands meet 
at the Fermi level. Close enough to Fermi level the band 
is nearly flat, which indicates the important role interac- 
tions might play in this material. 

C. Extracting hopping parameters from DFT 

Previously, bulk graphite (with the Bernal stacking 
order) SWM hopping parameters have been extensively 
studied using DFT and measured in experiments. How- 
ever, the values of the SWM parameters appropriate for 
ABC-stacked trilayer graphene were previously unknown. 
We extract their values by fitting the effective model with 
the DFT data in the zero electric field limit. The eigenen- 
ergies of the Hamiltonian in Eq. ^ in the absence of 
external potentials are 

E^^^ ^hs± ^Jhl^ + hi + 2cos(3v?p)/ich/itr , (7) 

where h^h = {voP)^/ji, Kr = 72/2 - 2uo«3P^/7i and 
hs = 5 — 2voVip^ /ji. To extract the remote hopping pa- 
rameters we first set the nearest neighbor in-plane hop- 
ping parameter 70 to 3.16 eV to set the overall energy 
scale. The values of S and. up to a sign, 72 can then be 
obtained by comparing the band energies at p = cal- 
culated by the two different methods. Then comparing 
£:(+) + ^(-) from the DFT data with Eq. (0), we obtain 
a value for 7470/71. Finally we notice that Eq. ([7]) im- 
plies that the gap between conduction (-|~) and valence 
(— ) bands vanishes at cos(3(y9p) = 1 if htr is negative and 
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FIG. 3: (Color online) The magenta curve is the DFT data 
while the Green one represents the effective model using the 
extracted parameters shown in Table|Il G — 47r/(\/3a) is the 
length of the reciprocal vectors and A; = is the K point. 

at cos(3(pp) = — 1 if htr is positive. Because of this prop- 
erty the Fermi level of a neutral balanced ABC trilayer 
is at the energy of three distinct Dirac points which are 
removed from the Dirac point separated in direction by 
27r/3. The triple Dirac point of the trilayer's simplified 
model is split into three separate single Dirac points. The 
DFT theory result that the conduction valence gap van- 
ishes along the K'M directions for which cos(3(^p) = 1 
implies that ht^ is negative and helps to fix the sign of 72. 
Values for 7370/71 and 7o/7i arc provided by the value 
of p at the Dirac points and the size of the splitting be- 
tween conduction and valence bands {2y/}i^^ + IiD along 
the cos(3v3p) = directions. The best overall fit we ob- 
tained to the bands around the K point and the deformed 
Dirac cones is summarized in Table HI where we com- 
pare with the corresponding fitting parameters for bulk 
graphitci^^ Our fit is extremely good in the low energy 
region in which we are interested, as shown in Fig. |31 
though there are still discrepancies as higher energies are 
approached. These discrepancies are expected because of 
the perturbative nature of the effective model and can be 
partly corrected by restoring the 1 — (woP)^/7i correction 
factor in Eq. ([S|). 



D. Electron(Hole) Pockets and Lifshitz Transitions 

With the effective model hopping parameters extracted 
from DFT we study the shape of the Fermi surface of a 
graphene trilayer. Fig. Hl^a) shows the constant energy 
contour plot of the electron band around zero energy. 
Clearly, under remote hopping the J=3 Dirac points 
evolve into three separate J = 1 Dirac points symmet- 
rically shifted away a little bit in the KM directions 
(/cx); each shifted Dirac point resembles a linear cone 
like the ones in monolayer graphene. The property that 
total chirality is conserved can be established by eval- 
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(a) Conduction band constant energy surfaces of an ABC 
graphene trilayer. 
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(b) Conduction band constant energy surfaces of an ABC 
graplicnc trilayer model with bulk graphite parameters. 
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FIG. 4: (Color online) Constant energy (in units of eV) con- 
tour plots of the conduction band near zero energy. G — 
47r/(\/3a) is the length of the reciprocal vector and = is a 
K point. (a)ABC-valued, (b)Bulk graphite valued Fermi sur- 
faces of a ABC trilayer. (a)The energies of the initial three 
electron pockets from inner to outer are 0.0, 2.5, 5.0, 6.0, and 
6.7 meV; The energies of the central triangles from outer to 
inner are 6.8, 6.9, 7.0, 7.1 and 7.2 meV; The energies of the 
bigger triangles from inner to outer are 6.8, 6.9, 7.0, 7.1, 7.2, 
7.5, 9.0, 10.0, 15.0, 20.0, and 30.0 meV. (b)The energies of 
the initial three electron pockets from inner to outer are 1.0, 
5.0, 7.5, 10.0, 10.2, 10.4 and 10.6 meV; The energies of the 
central triangles from inner to outer are 10.0, 10.2, 10.4 and 
10.6 meV; The energies of the bigger triangles from inner to 
outer are 10.8, 15.0, 20.0, and 30.0 meV. 



TABLE I: Summary of SWM hopping parameters obtained 
by fitting DFT bands in ABC-stacked trilayer graphene to a 
low-energy effective model. We compare with bulk graphite 
values from References.— 



parameters 


graphite(eV^) 


ABC trilayer (eV) 


5 


0.008 


-0.0014 


7i 


0.39 


0.502 


73 


0.315 


-0.377 


74 


-0.044 


-0.099 


72 


-0.020 


-0.0171 



uating Berry phases along circular paths far from the 
Dirac points where the remote hopping processes do not 
play an essential role. The Dirac point distortion occurs 
because the direct hopping 72 process does not involve 
2D translations and therefore gives a momentum inde- 
pendent contribution to the Hamiltonian which does not 
vanish at the Brillouin-zone corners. A similar distortion 
of the simplified-model ideal chirality Dirac point occurs 
in any 3m-layer system of ABC stacked (m is a positive 
integer) graphene sheets. Around each deformed Dirac 
cone there is a electron (hole)-like pocket in the conduc- 
tion(valence) band at low carrier densities and two Lif- 
shitz transitions^ as a function of carrier density. Take 
the conduction band for example. As shown in Fig.|4l^a), 
immediately above zero energy, the constant energy sur- 
face consists of three separate Dirac pockets. At the 
first critical energy 6.7 meV, the three electron pockets 
combine and a central triangle-like hole pocket appears. 
(Energies are measured from the Fermi energy of a neu- 
tral trilayer.) At this energy three band-structure saddle 
points occur midway between the shifted Dirac points, 
and thus the density-of-states diverges. Fermi levels close 
to these 2D logarithmic van Hove singularities could lead 
to broken symmetry states. At the second critical en- 
ergy 7.2 meV, the central pocket and the three remote 
pockets merge into a single pocket with a smoothed tri- 
angle shape. Fig. Sl^a) is in excellent agreement with 
constant energy surfaces constructed directly from our 
DFT calculations (Figure not shown). The two similar 
Lifshitz transition energies in the valence band occur at 
—7.9 meV and —9.9 meV. The constant energy surface 
at the second Lifshitz transition solves 

£;W(p^O) = i?(±)(p-0), (8) 

where +(— ) refers to conduction and valence band cases. 
This critical condition can be specified using the law of 
cosines as shown in Fig. [5l where for trilayers 0Borry = 
Stt and Hq = I72/2I ± 2woZ'4P^/7i. This momentum- 
dependent trigonometric condition can be easily gener- 
alized to the case of any other graphene multilayer and 
to the case with an external potential difference. Above 
the second Lifshitz transition, the constant energy sur- 
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FIG. 5: (Color online) A momentum-dependent trigonomet- 
ric relationship which describes how the shape of the con- 
stant energy surfaces near the Lifshitz transitions is collec- 
tively governed by chiral dispersion, trigonal warping, and 
Berry phases. 



face is triangular in shape, with a trigonal distortion that 
differs in orientation compared to the one obtained by 
plugging the bulk graphite values for the hopping pa- 
rameters into the same effective model Eq. © as illus- 
trated in Fig. m^b). The ABC-stacked trilayer trigonal 
distortion has a different orientation and is weaker. The 
difference mainly reflects a difference in the sign of 73, 
which favors anti-bonding orbitals at low energies. The 
warping of the constant energy surface becomes hexag- 
onal at 8 ^ 9 meV, which provides nearly parallel flat 
pieces on the edges of the hexagon leading to strong nest- 
ing. This might support some competing ground states 
and a density-wave ordered phase might then exist at 
a small but finite interaction strength. The electronic 
properties of low-carrier density systems in graphene tri- 
layers will be sensitive to these detailed band features. 
Future ARPES experiments should be able to determine 
whether or not these features arc predicted correctly by 
our DFT calculations. 



III. INDUCED BAND GAPS IN TRILAYERS 

A. Energy Bands with Electric Fields 

Fig. ini shows the energy band structure of a ABC- 
stacked graphene trilayer for several external electric po- 
tential differences between the outermost layers. In the 
presence of an external field, as in the graphene bilayer 
case, the energy gap is direct but, because the low-energy 
spectrum develops a Mexican hat structure as the electric 
potential difference increases, occurs away from the K or 
K' point. Charge transfer from the high-potential layer 
to the low-potential layer partially screens the external 
potential in both bilayer and multilayer cases. Fig. [7Ka) 
plots the screened potential U and Fig. [7l[b) the energy 
gap, as a function of the external potential Uext for bilay- 
ers and trilayers calculated using both DFT and the full 
band self-consistent Hartree approximation. The simple 
model Hartree calculations agree quite well with the DFT 
results generally. We find that the screening is stronger 
in a trilayer system, and that the maximum energy gap 
is slightly smaller. In both bilayer and trilayer, remote 



hopping suppresses the size of the energy gap but make 
little difference to the screening. 




k/G (direction: K ) 

FIG. 6: (Color online) The band structures of a ABC 
graphene trilayer with external electric potential differences 
between the outermost layers. The external potential dif- 
ference C/cxt values are O.O(red), 0.2(blue), l.O(green) and 
2.0(magenta) eV, respectively. G — 47v/{y/Sa) is the length 
of the reciprocal vectors and = is a K point. 



B. Self-Consistent Hartree Calculation 

As in the bilayer case, it is interesting to develop a 
theory of gap formation and external potential screening 
for ABC trilayers by combining the low-energy effective 
model with a Poisson equation which takes Hartree inter- 
actions into account. This simplified approach provides 
a basis for discussing the dependence on layer number for 
general N. We therefore consider an isolated graphene 
A'^-layer with an interlayer separation d = 0.335nm un- 
der an external electric field E^xt perpendicular to the 
layers, neglecting the finite thickness and crystalline in- 
homogeneity of the graphene layers. In an isolated sys- 
tem charge can only be transferred between layers so that 
n = rit + rib = 0. Defining Sn = Uh = —nt and using a 
Poisson equation, we find that the screened electric po- 
tential difference U between the outermost layers is 

U ^Ucxt + '^Tre^iN -l)d6n. (9) 

In the two-band effective model, 6n is accumulated 
through the layer pscudospin polarization of the valence 
band states and is thus given by the following integral 
over momentum space: 

7f^(^i(fc)l?|V^Kfe)), (10) 

Jbz l^""; ^ 

where the factor 2 accounts for spin degeneracy, |V'i(fc)) 
is a band eigenstate in the presence of E^xt, band in- 
dex i runs over all the filled valence band states, and 
az/'2 denotes the layer-pseudospin. Any Hamiltonian 
of a two-band model can be generally written as H = 
ho{p) + h{p) ■ (J. Defining tanflp = -^//i^ 4- h^/h-s and 
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(b)Energy gap evolution 
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FIG. 7: (Color online) Evolution of (a) the screened electric 
potential difTerence and (b) the energy gap, with respect to 
the increase of the external electric potential difference be- 
tween the outermost layers. □ represents the DFT calcula- 
tions while f (*) denotes the full band self-consistent Hartree 
calculations without (with) remote hopping 73. 



tan0p = h2/hi the conduction and valence band states 
in the sublattice representation are 




ill) 



dp I I cr^ 



27r fp, 



cos O-p p dp dif^ 



(12) 



where Pc — 71 /wo is the high momentum cutoff of the 
effective model and (pp is the angle of p. 

Let's first discuss the simplified two-band model which 
has only the chiral term. For general N 



H. 



(N) 
ch 



(-71)^-1 

(vopr 



(Trt)^ 




(-71)^ 



^(cos(7V(/7p)crx + sin(7V(/?p) cTy)(13) 



The electric potential in the two-band model is ±-^^ ctz 



Inserting Eq. in Eq. ^ and Eq. we obtain an 
algebraic formula for the self-consistent Hartree potential 
valid for general N: 



C/oxt U , 4{N-l)dm2 



71 



F{N,U) = - 



71 

1 '•'^ 



ao 
dt 



FiN, U) , (14) 



2^,(1 i -(^)^) 

2^ IV ATI r,' AT ' ^ TT ) ) ■ 



^iV'2' N 



U 



where ag = 0.053nm is the Bohr radius, TO2 is the effec- 
tive mass of a graphene bilayer, tc = (2ji/U)'^/^ and 2F1 
is Gauss' hypergeometric function. In the limit of large 
N, F{N, J7) — )- 1 and thus the Hartree equation reduces 
to 



4{N -l)dm2 

U ~ [/oxt 71 



(15) 



except at very small U . For small U and N ^ 2, the 
Hartree equation reads 



f/c: 



u 



2d 1712 471 
ao mo U 



(16) 



which is consistent with previous Hartree calculations in 
graphene bilayersi^. In the limit of small J7 for TV > 2, 
the Hartree equation has the asymptotic form 



U 



271 V 271 

where the factor C 



C, 



(17) 

1 ^)]-^/2, 



2(Af-l)d m2 (-1 2 

ao m„ 2-N~' 2-3N ' 

The larger the value of N, the flatter the chiral bands, 
and the stronger the screening. For N = 2 the screen- 
ing response is linear up to a logarithmic factor, while for 
larger N, superlinear screening leads to a screened poten- 
tial difference which initially grows slowly with external 
potential following U oc U^(^ . The strongest possible 
screening reduction of the external potential corresponds 
to the Hartree-potential due to transfer of all the states 
in the energy regime < 271 over which the low energy 
model applies to one layer. 

For the trilayer case we can perform a similar calcu- 
lation using the full low-energy Hamiltonian derived in 
Eq. ([6]). In this case we find that 



Uext 


U 


71 


71 


G{U) = 





8d 7712 

ao me 



(18) 



dt- 



(h,h + ^/2ht,)^ + hl 



-) 



t^,h,, = \^-'-^tlh,,p = ^il-t), 
t = {^y^, and K{x) is the complete elliptic integral of 
the first kind. Fig. [5] compares the screening properties of 



where /ich 

'VoP\2 
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the full low-energy effective model for trilayers to the chi- 
ral model results for TV = 2, 3, 4, 5. For ?7ext < 7i/2, the 
energy regime over which the low-energy effective model 
applies, we see that screening increases systematically 
with N because of smaller gaps between conduction and 
valence band orbitals which make the occupied valence 
band orbitals more polarizable. The comparison between 
the simplified chiral model and the low-energy effective 
model for = 3 demonstrates that remote hopping pro- 
cesses suppress screening because they tend to increase 
the gap between conduction and valence bands at mo- 
menta near the Brillouin-zone corner. 
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FIG. 8: (Color online) U v.s. (7cxt plot describes the screening 
effect in different chiral- A*" systems. The Chiral model results 
refer to the Hamiltonian in Eq. (|13p while the full model re- 
sults refer to the Hamiltonian in Eq. 

In concluding this section we caution that occupied 
a orbitals, neglected in the low-energy effective model 
and TT-band tight-binding models, will contribute slightly 
to polarization by an external electric field and there- 
fore to screening. Furthermore exchange potentials will 
also be altered by an external electric field and infiuence 
the screening. Since exchange interactions are attrac- 
tive, they always work against screening and will make 
a negative contribution to the screening ratio we have 
discussed in multilayers. Because the low energy eigen- 
states in multilayers are coherent superpositions of states 
localized in different layers, our DFT calculations which 
employ a local exchange approximation, might also yield 
inaccurate results for the screening ratio. In fact simple 
measurements of the screening properties might provide a 
valuable window on many-body physics in ABC-stacked 
graphcne multilayers which lies outside the scope of com- 
monly employed approximations. 

IV. DISCUSSION 

We have derived an effective model for the low- 
energy conduction and valence bands of an ABC-stacked 
graphene multilayer. The low-energy model can be 
viewed as a momentum-dependent pseudospin Hamilto- 



nian, with the pseudospin constructed from the low en- 
ergy sites on the top and bottom layers. The simplified 
version of this model starts from a 7r-band tight-binding 
model with only nearest neighbor hopping and yields a 
pseudospin magnetic field whose magnitude varies as mo- 
mentum in an N-layer stack and whose direction is 
A^^p where (f>p is the momentum orientation. The likely 
importance of electron-electron interactions in multilay- 
ers can be judged by comparing the characteristic band 
and interaction energies in a system with carrier den- 
sity n and Fermi wavevector pp cc \/n. The character- 
istic Coulomb interaction energy per-particle in all cases 
goes like e^7i^/^, while the band energy goes like n^/^. 
For low-carrier densities the band energy scale is always 
smaller. In the case of trilayer ABC graphene, the inter- 
action energy scale is larger than the band energy scale 
for carrier density n < 10^'^cm^^. 

Although interactions are clearly important and can 
potentially introduce new physics, the chiral band model 
is not valid at low-densities because of the influence of 
remote hopping processes which we have estimated in this 
article by carefully fitting a low-energy effective model 
to DFT bands. The Hamiltonian in Eq. © combined 
with the parameters in Table I should be used to describe 
graphene trilayers with low carrier densities. In a realistic 
system the Fermi surface of a ABC trilayer with a low 
carrier density consists of three electron pockets centered 
away from the K point. As the carrier density grows 
these pockets convert via a sequence of two closely spaced 
Lifshitz transitions into a single K-ccntercd pocket. The 
carrier density at the Lifshitz transition is ~ lO^^cm^^, 
which translates to a Coulomb interaction scale of ^ 45 
meV, compared to a Fermi energy of ~ 7 meV. 

The Berry phase associated with the momentum- 
dependence of the pseudospin orientation field, tt for a 
full rotation in single-layers and 27r in the bilayer chiral 
model for example, is known^"— to have an important 
infiuence on quantum corrections to transport. Because 
of their very different Berry phases time-reversed paths 
are expected to interfere destructively for iV-odd systems 
while constructively for TV-even system, leading to weak 
anti-localization for odd TV and weak localization for even 
TV. This general tendency will however be altered by 
trigonal and other corrections to the low-energy effective 
Hamiltonian, like those we have derived for trilayers. The 
infiuence of these band features on quantum corrections 
to transport can be evaluated starting from the results 
obtained here. 

Another important consequence of Berry phases in the 
chiral model is the unconventional Landau level struc- 
ture it yields^i^"— . In the chiral model for ABC trilay- 
ers there is a three- fold degeneracy at the Dirac point, 
in addition to the usual spin and valley degeneracies. 
This grouping of Landau level leads to the expecta- 
tion that quantum Hall studies in trilayers will reveal 
plateaus that jump from one at —6e^/h to one at 6e^ /h. 
Electron-electron interactions acting alone are expected 
to lift these degeneracies and give rise to quantum Hall 



ferromagnetismiSri^. These interaction effects will act in 
concert with small corrections to the Landau level struc- 
tures due to the remote hopping terms that have been 
quantified in this paper. 

Although wc have discussed the case of ABC stacked 
trilayers, we expect qualitatively similar results for ABC 
stacking sequences of general thickness TV. At low ener- 
gies the band structure will consist of a conduction and 
a valence band with dispersion and a gap in the pres- 
ence of an external electric field across the film. In the 
presence of a magnetic field N Landau levels are pinned 
to the neutral system Fermi level for each spin and valley. 
At the lowest energies, within around lOmeV of the neu- 
tral system Fermi level, constant energy surfaces will be 
strongly influenced by remote hopping processes which 
will also split the Dirac point Landau levels. The re- 
mote hopping terms give rise to saddle-points in the band 
structure at which the density-of-states will diverge. Bro- 
ken symmetry electronic states arc mostly likely to oc- 
cur when the Fermi level is coincident with these saddle 
points. The energy range over which the low-energy ef- 
fective model applies will, however, decrease with film 
thickness. We expect both disorder and interaction ef- 
fects to be strong within this family of low-dimensional 
electron systems, which should be accessible to experi- 
mental study in samples for which disorder is weak on the 
energy scale over which the low-energy effective model 
applies. 

In summary, we have derived an effective model for 
trilayers, extracted the hopping parameters for ABC- 
stacked multilayers, from DFT and studied the trilaycr 
Fermi surfaces. Furthermore, we have explored the 
screening effect in trilayers and then explained and com- 
pared with other C2DES cases by a tight-binding model 
self-consistent Hartree method. Lastly, we have ar- 
gued the importance of Berry phases and interactions 
in C2DES. 
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FIG. 9: (Color online) Schematic of hoppings from Ai to Ba; 
(a) one-step Ai — >■ Ba and (b) three-step Ai — >■ B1A2 — >■ 
B2A3 — >• B3 and (c) (d) two-step Ai — >■ B1A2 — > B3 and 
Ai — > B2A3 — > B3. Schematic of hoppings from Ai to Ai; (e) 
two-step Ai — > -B1A2 — > A\ and (f) two-step A\ — i> A2B\ — >• 
^1. 



Appendix A: Diagrammatic Derivation of the Low 
Energy Effective Model 
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As a result of tight-binding model, each term of the 
effective Hamiltonian Eq. (|6]) has a unique physical pic- 
ture. Hereafter, we view the strongly stacked pair BiAi+i 
as a single dimer site and assume zero external poten- 
tials for simplicity. The general formula of effective low 
energy models Eq. ([5]) can be understood as following. 
The terms in the second parenthesis represent the lead- 
ing hopping processes, while the terms in the first paren- 
thesis are approximately 1 — (wop/71)^ and give a small 
correction. Hn is the unperturbed Hamiltonian of low 
energy sites and thus includes the direct hopping and on- 
site energy. H21 and H12 are hoppings from and to low 
energy sites, respectively, describing the coupling to high 
energy ones. H22 contains the hoppings between high 
energy sites and is an intermediate process. Therefore 
Hi2{H22)~^ H2\ together gives the general "three"-step 
hoppings which start from and end at low energy sites 



10 



by way of high energy ones. Note that the intermediate 
process within high energy sites is zero for single lay- 
ers, a constant for bilayers, one-step for trilayers, and 
multi-step for > 3 layers. In bilayers for example, the 
linear trigonal warping term arises from Hn, while the 
chiral term attributes to Hi2{H22)^^ H2i- Because H22 
gives no hopping and is simply 71, Hi2{H22)~^ H21 is re- 
duced to two-step and hence the chiral term is quadratic. 
In the trilayer case, for the matrix element B^Ai, Hn 
provides the first term of i?tr shown in Fig. [HJa) while 
Hi2{H22)~^ H21 contributes Hch and the second term 
of iJtr as depicted in Fig. IHl^b) and (c) (d) , respectively. 
Hi2{H22)~^ H21 also gives rise to the second term of Hg 
for the matrix element AiAi as presented in Fig. |9l[e)(f). 
Generally, in order to derive the low energy effective 



model for a general ABC-stacked TV-layer graphene, we 
first need to write a 27V x 2N Hamiltonian matrix as 
Eq. (121) , then we specify all the leading hopping processes 
in the diagrammatic language like Fig. 1^1 instead of in- 
verting the large Hamiltonian matrix. The hopping di- 
agrams are convenient for systematic calculations in a 
way similar to the way Feynman diagrams help in per- 
turbation theories. The exact coefficient of one hopping 
process can be easily calculated using Eq. ([5|) by picking 
up the starting and ending sites, setting matrix elements 
of unrelated sites as zero and turning off the unrelated 
hopping parameters. Frequently, one hopping process 
can be neglected because its requirement of more than 
one sub-hopping with comparably small amplitudes. 
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